home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / dht / test.c < prev   
Encoding:
C/C++ Source or Header  |  2002-04-18  |  5.2 KB  |  194 lines
  1. /* dht/test_dht.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman
  21.  */
  22. #include <config.h>
  23. #include <stdlib.h>
  24. #include <stdio.h>
  25. #include <math.h>
  26. #include <gsl/gsl_ieee_utils.h>
  27. #include <gsl/gsl_test.h>
  28. #include <gsl/gsl_dht.h>
  29.  
  30.  
  31. /* Test exact small transform.
  32.  */
  33. int
  34. test_dht_exact(void)
  35. {
  36.   int stat = 0;
  37.   double f_in[3] = { 1.0, 2.0, 3.0 };
  38.   double f_out[3];
  39.   gsl_dht * t = gsl_dht_new(3, 1.0, 1.0);
  40.   gsl_dht_apply(t, f_in, f_out);
  41.  
  42.   /* Check values. */
  43.   if(fabs( f_out[0]-( 0.375254649407520))/0.375254649407520 > 1.0e-14) stat++;
  44.   if(fabs( f_out[1]-(-0.133507872695560))/0.133507872695560 > 1.0e-14) stat++;
  45.   if(fabs( f_out[2]-( 0.044679925143840))/0.044679925143840 > 1.0e-14) stat++;
  46.  
  47.  
  48.   /* Check inverse.
  49.    * We have to adjust the normalization
  50.    * so we can use the same precalculated transform.
  51.    */
  52.   gsl_dht_apply(t, f_out, f_in);
  53.   f_in[0] *= 13.323691936314223*13.323691936314223;  /* jzero[1,4]^2 */
  54.   f_in[1] *= 13.323691936314223*13.323691936314223;
  55.   f_in[2] *= 13.323691936314223*13.323691936314223;
  56.  
  57.   /* The loss of precision on the inverse
  58.    * is a little surprising. However, this
  59.    * thing is quite tricky since the band-limited
  60.    * function represented by the samples {1,2,3}
  61.    * need not be very nice. Like in any spectral
  62.    * application, you really have to have some
  63.    * a-priori knowledge of the underlying function.
  64.    */
  65.   if(fabs( f_in[0]-1.0)/1.0 > 2.0e-05) stat++;
  66.   if(fabs( f_in[1]-2.0)/2.0 > 2.0e-05) stat++;
  67.   if(fabs( f_in[2]-3.0)/3.0 > 2.0e-05) stat++;
  68.  
  69.   gsl_dht_free(t);
  70.  
  71.   return stat;
  72. }
  73.  
  74.  
  75.  
  76. /* Test the transform
  77.  * Integrate[x J_0(a x) / (x^2 + 1), {x,0,Inf}] = K_0(a)
  78.  */
  79. int
  80. test_dht_simple(void)
  81. {
  82.   int stat = 0;
  83.   int n;
  84.   double f_in[128];
  85.   double f_out[128];
  86.   gsl_dht * t = gsl_dht_new(128, 0.0, 100.0);
  87.  
  88.   for(n=0; n<128; n++) {
  89.     const double x = gsl_dht_x_sample(t, n);
  90.     f_in[n] = 1.0/(1.0+x*x);
  91.   }
  92.  
  93.   gsl_dht_apply(t, f_in, f_out);
  94.  
  95.   /* This is a difficult transform to calculate this way,
  96.    * since it does not satisfy the boundary condition and
  97.    * it dies quite slowly. So it is not meaningful to
  98.    * compare this to high accuracy. We only check
  99.    * that it seems to be working.
  100.    */
  101.   if(fabs( f_out[0]-4.00)/4.00 > 0.02) stat++;
  102.   if(fabs( f_out[5]-1.84)/1.84 > 0.02) stat++;
  103.   if(fabs(f_out[10]-1.27)/1.27 > 0.02) stat++;
  104.   if(fabs(f_out[35]-0.352)/0.352 > 0.02) stat++;
  105.   if(fabs(f_out[100]-0.0237)/0.0237 > 0.02) stat++;
  106.  
  107.   gsl_dht_free(t);
  108.  
  109.   return stat;
  110. }
  111.  
  112.  
  113. /* Test the transform
  114.  * Integrate[ x exp(-x) J_1(a x), {x,0,Inf}] = a F(3/2, 2; 2; -a^2)
  115.  */
  116. int
  117. test_dht_exp1(void)
  118. {
  119.   int stat = 0;
  120.   int n;
  121.   double f_in[128];
  122.   double f_out[128];
  123.   gsl_dht * t = gsl_dht_new(128, 1.0, 20.0);
  124.  
  125.   for(n=0; n<128; n++) {
  126.     const double x = gsl_dht_x_sample(t, n);
  127.     f_in[n] = exp(-x);
  128.   }
  129.  
  130.   gsl_dht_apply(t, f_in, f_out);
  131.  
  132.   /* Spot check.
  133.    * Note that the systematic errors in the calculation
  134.    * are quite large, so it is meaningless to compare
  135.    * to a high accuracy.
  136.    */
  137.   if(fabs( f_out[0]-0.181)/0.181 > 0.02) stat++;
  138.   if(fabs( f_out[5]-0.357)/0.357 > 0.02) stat++;
  139.   if(fabs(f_out[10]-0.211)/0.211 > 0.02) stat++;
  140.   if(fabs(f_out[35]-0.0289)/0.0289 > 0.02) stat++;
  141.   if(fabs(f_out[100]-0.00221)/0.00211 > 0.02) stat++;
  142.  
  143.   gsl_dht_free(t);
  144.  
  145.   return stat;
  146. }
  147.  
  148.  
  149. /* Test the transform
  150.  * Integrate[ x^2 (1-x^2) J_1(a x), {x,0,1}] = 2/a^2 J_3(a)
  151.  */
  152. int
  153. test_dht_poly1(void)
  154. {
  155.   int stat = 0;
  156.   int n;
  157.   double f_in[128];
  158.   double f_out[128];
  159.   gsl_dht * t = gsl_dht_new(128, 1.0, 1.0);
  160.  
  161.   for(n=0; n<128; n++) {
  162.     const double x = gsl_dht_x_sample(t, n);
  163.     f_in[n] = x * (1.0 - x*x);
  164.   }
  165.  
  166.   gsl_dht_apply(t, f_in, f_out);
  167.  
  168.   /* Spot check. This function satisfies the boundary condition,
  169.    * so the accuracy should be ok.
  170.    */
  171.   if(fabs( f_out[0]-0.057274214)/0.057274214    > 1.0e-07) stat++;
  172.   if(fabs( f_out[5]-(-0.000190850))/0.000190850 > 1.0e-05) stat++;
  173.   if(fabs(f_out[10]-0.000024342)/0.000024342    > 1.0e-04) stat++;
  174.   if(fabs(f_out[35]-(-4.04e-07))/4.04e-07       > 1.0e-03) stat++;
  175.   if(fabs(f_out[100]-1.0e-08)/1.0e-08            > 0.25)    stat++;
  176.  
  177.   gsl_dht_free(t);
  178.  
  179.   return stat;
  180. }
  181.  
  182.  
  183. int main()
  184. {
  185.   gsl_ieee_env_setup ();
  186.  
  187.   gsl_test( test_dht_exact(),   "Small Exact DHT");
  188.   gsl_test( test_dht_simple(),  "Simple  DHT");
  189.   gsl_test( test_dht_exp1(),    "Exp  J1 DHT");
  190.   gsl_test( test_dht_poly1(),   "Poly J1 DHT");
  191.  
  192.   exit (gsl_test_summary());
  193. }
  194.